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I. INTRODUCTION 


The extraction of ocean tidal components from the analysis of satellite altimetry data could become 
an alternative and complementary method to the numerical integration of the Laplace tidal equations. 

Mazzega (1985) created a global model of the M2 tide using twenty-four days of SEASAT altimetry 
data; the solution was obtained by means of surface spherical harmonics; the results were qualitatively re- 
alistic. No hydrodynamic equations were used. 

Woodworth, P.L. and Cartwright, D.E. (1986) have extracted the M2 ocean tide from SEASAT al- 
timetry data. They used three complementary methods. The first method provides point measurements of 
the tide at the crossovers of the SEASAT repeat orbit ground track; it was applied in the tropical ocean 
areas. The other two methods involve the spatial expansion of M2 in terms of either surface spherical har- 
monics or Platzman normal modes of the world ocean. The results reproduce many features of the tide 
represented in recent tidal models. 

The purpose of this investigation is to develop an estimation technique which will serve to extrap- 
olate tidal amplitudes and phases over entire ocean basins using existing gauge data and the precise alti- 
metric measurements which are now beginning to be provided by satellite oceanography. The applicability 
of the extrapolating technique was tested in the Lake Superior basin by Sanchez, Rao and Wolfson (1985) 
and in the Atlantic-Indian Oceans(using a 6° x 6° degree) by Sanchez, Rao and Steenrod (1986 a, b). 

The method to be used in this investigation requires several distinct steps. First it is necessary to de- 
termine numerically the stream function and velocity potential orthogonal functions (The Stokes/Helm- 
holtz Potentials) which span the space of the basin under consideration. These space functions are then 
used in the Laplace tidal equations to determine the homogeneous solution (normal modes) and the 
forced solution. The latter is obtained by adding the astronomical forcing function modified to include 
solid-earth tides. 

The velocity potential eigenfunctions obtained as a first step are used also to extrapolate the surface 
height field over the entire space domain of the given basin and this approach will constitute a distinct 
and integral part of the investigation. 

The theoretical foundation is Proudman’s theory (1918) as formulated by Rao (1966). The theory 
provides the formalism for calculation of the gravitational (first class) normal modes and the rotational 
(second class or Rossby waves) normal models of irregularly shaped basins with realistic bathymetry. 

The method requires the solution of two elliptic partial differential equations with second order oper- 
ators which are simpler than the tidal operator. The boundary conditions correspond to vanishing of the 
stream function and normal derivative of the velocity potential. The elliptic operators are represented nu- 
merically in finite difference form, the grid used is a Richardson lattice which preserves self-adjoint ness. 
The solutions yield the velocity and surface height fields in terms of orthogonal functions with time-depen- 
dent coefficients. These functions are then substituted into the Laplace’s tidal equations: if the homo- 
geneous equations are used one obtains the normal modes; if the forcing terms are included then the 
forced solution is obtained. In both cases the solution is obtained numerically. The surface height field is 
only dependent on the velocity potential orthogonal functions. The expansion coefficients of these func- 
tions can be estimated in a least-square sense from available selected tidal measurements. 
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2. BASIC EQUATIONS 


Free Solutions 

The method of approach was developed originally by Proudman (1918) using a Lagrangian approach. It was reformulated 
by Rao (1966) from the Eulerian point of view. The basic ideas of the method presented below follow Rao’s line of devel- 
opment. The basic equations are the linearized shallow water equations on a rotating plane: 


♦ 

™ - flM] = -g Hh Vr, 

4^ + v • M = 0 
at 

where 


( 2 . 1 ) 


M = HV = (M,N) 


V s (u,v) 
f = 2 co sin 0 
h(x,y) s H (x,y)/H. 

■+ 

H(x,y) is the variable depth of the fluid in equilibrium and H is some constant scaling depth, f is the Coriolis parameter, V is 
the horizontal velocity vector, r\ the fluctuation of the free surface, g is the apparent gravitational acceleration and V is the 
horizontal gradient operator. [] : denotes rotation of the vector through a right angle in the clockwise direction of the hori- 
zontal plane, i.e. [] = — (k x V), k being a unit vector vertical to the horizontal plane. 

The appropriate boundary conditions to be adjoined to equations (2.1) are 

M * n = 0 (2.2) 

-*■ 

on the coast, where n is the unit normal to the coastline. 

The transport vector M may be partitioned as follows: 

M = + M v (2.3) 

where 

M* = -h V4> 

iCr = -[v«t>] (2.4) 


M 4 ' is the solenoidal part of M while h is the irrotational part since 

V'[h -1 M^] = 0 


vr = o 

To complete the determination of >}> and ip, it is necessary to specify the boundary conditions ' n = 0 and M v • n = 0 to 
ensure that Eq. (2.2) is satisfied. In terms of <(> and v the conditions then are 

h = 0 (2.5) 

an 
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and 


ip = 0 


on the boundaries. The divergence of the transport field and the vorticity of the velocity field yield 

V • M = -VhV* 

V • [h -1 M] = Vh~ l Vip (2.6) 

-► 

If M is known as a function of the horizontal coordinates, the left-sides of (2.6) are specified. Then each equation represents 
an inhomogeneous elliptic differential equation with homogeneous boundary conditions given by Eq. (2.5), and^it is well 
known that such problems possess a unique solution. It is also straightforward to prove that the representation of M as given 
in Eqs. (2.3, 4, 5) is unique. Since M itself is unknown, the procedure then consists of expanding * and ip in terms of the 
spectra of the elliptic operators appearing in Eqs. (2.6); that is, we seek solutions which satisfy the following equations 


V • hV<j> y = -A y * 

V • h~'Vl^ r = -fiylf/y 


(2.7) 


where * y , ip Y are the characteristic functions and A y , p Y are the characteristic values associated with the corresponding op- 
erators Vh V and V h" l V. The characteristic functions satisfy the boundary conditions 



0 


h l \Py = 0 


(2.8) 


The condition h~“*ip y = 0 imposes a stronger condition than that required by (2.5). However, the factor h“* is necessary to 
make the ip Y problem self-adjoint. Since the problems (2.7, 8) are self-adjoint, the characteristic values of A y and p Y are real 
and the characteristic functions * y and \p Y each constitute a complete and internally orthogonal set. Without loss of gener- 
ality, the orthogonality condition may be chosen as: 


/h _1 M+-M^dA = A y J<t» y tyl A = Ac 2 H 2 d yP , 

= XyfyyittpdA = Ac 2 H 2 d yP , 


(2.9) 


where c 2 = gH and A is the surface area of the basin. 6 y p is the Kronecker delta. We have further defined, in accordance with 
(2.4), 


M+ = -hV$ y and M* = -[V Vy ] 


in the orthogonality condition (2.9). 


The components of the transport vector can now be represented by the sums 

M* = I P y M* 

M* = I Q r M£ 


( 2 . 10 ) 
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where P y and Q y are the expansion coefficients. The orthogonality conditions ensure that equations (2.10) represent the 
least squares approximations to and if and when the sums span the complete spectra of equations (2.7). The height 
field r 7 is governed by the divergent part of M and the $ y functions yield a sufficient basis for its representation. A conveni- 
ent representation for t? may be taken as: 


n = P/v 
n r = c- l a Y )'% , 


(2.11) 


where R y are the expansion coefficients for the 17 -field. 

Substitution of equations (2.10) and (2.11) into equations (2.1) and the use of the orthogonality conditions (2.9) yield a set 
of prediction equations for the expansion coefficients: 


—y-1 AypPp- 1 BypQp-OyRy = 0 
dt ft p 

7 , -fVr| D ,iQr» 

dt P P 


where the coupling coefficients are given by 

A yp = {My, [Mj]} , 

Cyp = {M}f, [Mj]} , 


B y/3 = {M}. [MJ]} , 
D y p = {M*. [MjJ ]} , 


(2.13) 


The quantity inside the angular brackets represents an inner product of two vectors W and T and is denoted by 

{W.T} = (1/c 2 AH 2 ) / A fh _1 WTdA. 


It can be seen from Eq. (2.12, 13) and the definition of the inner product that all the coupling coefficients vanish when the 
Coriolis parameter f = 0 and hence 


+ *2 rPy = o. 

dt 2 

i l y = (c 2 ky ) 1/2 , is the non-rotating frequency of oscillation. In equations ( 2 . 12 , 13) the subscript y (or p) is used as a proxy 
for a binary index and represents an ordering of the characteristic functions <j> y , ip y in some as yet unspecified manner. For 
convenience we replace the wavenumber index 7 by scalar indices i = 1 ,2,3 (or j = 1 ,2,3 . . .) and denote 

Aj = Ayj? — Myj, ilj — $ y i 
Pj = Pyj, = Qyj, Rj = Ry| 

A ij = Ayi^j 

Bij = B ri pj (2.14) 

Cij = 

°ij 3 D r i/3j 
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Using the above notation, we can define column vectors 


P = col (Pi), Q = col (Qj), R = col (Rj) 


and matrices 


S ■ 


+ 

P 

Q 


R 


A S { Ay }, B S { Bij } 
£ s { Cjj }, D s { Dy } 


<$> = diagonal i?i 


Equations (2.12) may now be written in the form 


where a is the square matrix 

no 


a 

no 


• • 


+ 

a s = 

0 

no 

-A 

-B 

-<£> 

-C 

-D 

0 

<0> 

0 

0 


In seeking the solutions for the normal modes, we assume 

s ~ e iot 

where a is the normal mode frequency with rotation and i = ✓ - 1 . Eq. (2.15) then reduces to 

(oil - ia ) s =0. 


(2.15a) 


(2.15b) 


(2.16) 


(2.17) 


(2.18) 


In the preceding equation II is the identity matrix, o’s are the characteristic values of the matrix ia. From the definition 
of the coupling coefficients given in (2.13) and the matrix a as defined in (2.17), it is clear that ia^ exhibits Hermitian sym- 
metry and hence the o’s constitute a spectrum of real eigenvalues. In computing the matrix elements in a^ , the basis func- 
tions 4>i and ipj are chosen to correspond to an ordering of the characteristic values Aj and \i x arranged so that Aj < A2 < A3... 
and fji\ < \i 2 < M3 * Such an ordering has been chosen since the Aj’s and ^’s have the dimensions of (wavenumber) 2 . Hence, 
at any order of truncation, those <j>j’s and tpj’s with the largest space scales are taken into account. 

Forced Solution 

In the computation of the forced solution we include the effects of the yielding of the solid earth to tide generating forces. 
The effects due to self-attraction of the tide and tidal loading as well as the frictional effects have not been included. The 
theory allows for the inclusion of these effects but they introduce computational complications and physical uncertainties 
and they are not necessary in the context of the application of the interpolation technique. Equations (2.1) are recast for 
the forced tidal oscillations as 
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(2.19) 



3. ESTIMATION 

The orthogonal functions <|> y form a basis for the expansion of the height field 17, as discussed earlier. Since these functions 
are characteristic for a particular basin and are described at all points inside the basin (within the resolution of the finite dif- 
ference grid), they can be used as optimal functions for extrapolation of the tidal field over the whole basin, given data at 
some selected points. 

Eventhough the basis functions are orthogonal, in fitting the data to these functions it is preferable to use linear least- 
squares techniques to determine the coefficients of expansion rather than the orthogonality property since real data are not 
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usually available at regularly spaced points. A similar approach was taken in a previous investigation by Rao and Schwab 
(1981) in which they determined the steady circulation in a closed basin for which the appropriate orthogonal functions are 
the xp y functions and Sanchez et al. (1985) used the <j> y functions for tidal extrapolation in Lake Superior. An outline of the 
procedure and the basic equations involved are given below. 

The tidal height field as given by equations (2.11) and (2.15) can be written 

r](x,y,t) = IR yr rj y (x,y) cosdt - IR yI rj Y (x>y) sindt (3.1) 

where R yr and R y j denote the real and imaginary parts of R y . The tidal height field can be expressed also in terms of ampli- 
tude A(x,y) and phase 0(x,y), that is 


ry(x,y,t) = A(x,y) cos [ot - 0(x,y)] 


(3.2) 


Comparing equations (3.1) and (3.2) one obtains the following relations, 


A(x,y) = [(I R yr n r ) 2 + (^ R yl n Y ) 2 ] U2 
0(x,y) = arc tan [-(ZR yI yj y )/(ZR yr yj y )] 


(3.3) 


Also, from equations (3.1) and (3.2) one obtains the following vector-matrix relations 


[rj ] R r = A cos 0 

~ # ^ (3.4) 

[?7 ] Rl = A sin 0 

The components of the vectors A cos 0 and A sin 0 are the “n” available tidal measurements. The matrix [i?] is available 
from the solution of the velocity potential eigenfunctions, its “k” columns will correspond to the velocity potential eigen- 
functions chosen to represent the forced solution, its “n” rows correspond to the location of the tidal measurements in the 
basin under consideration. The least squares solution to equations (3.4) is then given by 


Rr = (fo] T [n]r 1 fo] T A COS 0, 

<\j f\J f\l 

Ri = ([nl T h]) _1 h] T A sin 0, 

r\j r\j r\j 


(3.5) 


rp 4 ▼ 

where [ry] 1 is the transposed matrix. Having determined a certain number of coefficients R r and Rj, since the functions 
ry y (x,y) are known over the entire basin, one can obtain the amplitude and phase of the tidal height at all points in the basin 
from eq. (3.3). The maximum number of coefficients fe r and fej that can be determined will be the same as the number of 
observations available. 


4. RESULTS 


Normal Modes and Forced Solutions 

A normal mode solution for the Atlantic and Indian Oceans has been obtained by means of a 3° x 3° finite difference grid. 
There are 1915 velocity potential and 1717 stream function points distributed such as to form a single Richardson lattice. 
A finite difference solution of equations (2.7) subject to the boundary conditions given by equations (2.8) yields the eigen- 
values and eigenvectors for both fields. The eigenfunctions of the 0 solution represent the non-rotating gravitational normal 
modes of the basin, the periods of oscillation of the lowest modes are given in the first column of Table 1 . The stream function 
modes generate the vorticity component of the flow field in the rotating case for the gravitational and rotational species, they 
are more dominant for the quasi-geostrophic rotational modes than for the gravitational modes. The normal mode solution 
was obtained by including the lowest 150 eigenfunctions from each field (0 and 0) into the dynamical equations, the solution 
of the eigenvalue problem posed by equation (2.18) yields the normal modes. The normal modes fall into two distinct cate- 
gories: the inertia -gravitational modified by rotation and the rotational modified by divergence. If the secular determinant is 
truncated at a size n x n, then one obtains 2n/3 gravitational modes and n/3 rotational modes, in pairs (plus and minus the 
same value). In our case there are 300 gravitational modes and 150 rotational modes. The rotational modes are characterized 
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by high ratios of kinetic to potential energy; these modes will vanish in a basin of constant depth and Coriolis parameter. The 
rotational normal modes are extremely sensitive to the resolution of the shape of the basin and the bottom topography and 
the convergence of their roots is more complicated than for gravitational modes. A detailed discussion of the nature of the 
two types of modes is given by Platzman (1975) and Rao and Schwab (1976). The gravitational modes are by far the most 
important in the context of the diurnal and semidiurnal tides as shown by Platzman (1984). 

Table 1 presents a list of the first 20 gravitational modes in order of decreasing period. Both the non-rotational and ro- 
tational values are given, the non-rotational gravity modes are equivalent to the eigenfunctions of the velocity potential or 




TABLE 1. GRAVITATIONAL MODES 


PERIOD (hours) 
NON-ROTATIONAL 

PERIOD (hours) 
ROTATIONAL 

ENERGY RATIO 
(KINETIC/POTENTIAL) 

3x3 

6° x 6° 

">0 I ® 

3x3 

6° x 6° 

3x3 

6 x6 

81.20 

77.68 

68.14 

67.92 

1.66 

1.47 

42.19 

40.98 

43.99 

42.53 

0.93 

0.96 

29.23 


31.16 


0.79 


29.03 

29.00 

29.20 

29.61 

0.94 

0.81 

28.21 

22.17 

28.20 

23.55 

0.94 

0.84 

22.13 


22.89 


0.90 


20.79 

20.73 

20.61 

20.93 

0.87 

0.86 

19.31 

19.20 

19.13 

19.40 

1.20 

1.19 

18.94 


17.59 


0.99 


16.29 

16.52 

16.21 

17.34 

0.92 

1.03 

15.34 


15.51 


1.05 


14.76 

14.56 

14.79 

15.73 

0.89 

1.11 


14.17 


14.44 


0.86 

14.05 

13.42 

13.68 

13.93 

0.98 

1.03 

13.13 

12.94 

12.85 

12.99 

1.05 

1.07 

12.89 

12.74 

12.68 

12.55 

1.05 

1.27 

12.25 


12.05 


1.15 


12.18 


11.88 


1.03 


11.73 

11.72 

11.55 

11.99 

1.03 

1.09 

11.57 

11.47 

11.10 

11.13 

1.01 

1.14 

11.30 


10.96 


1.11 







Proudman functions. The last entry is the ratio of energies (kinetic/potential). Table 1 gives the values obtained with the 
3° x 3° grid as well as the results previously obtained with a 6° x 6° grid (1986 a, b). Notice that in the range covered there 
are only 14 values for the 6° x 6° grid as compared to 20 for the finer 3° x 3° mesh; this is to be expected. 

Table 2 lists the 20 most powerful modes for the M2 and K1 tidal components in order of decreasing power. The entries 
in the table give the power ranking, the rotational period in hours, the percentage of the total power for that mode and the 
percentage of the total power contributed by all the modes included up to that point. 


TABLE 2. POWER SPECTRUM FOR THE M2 AND K1 COMPONENTS. 
M2 K1 


RANK 

PERIOD 

(HOURS) 

POWER (%) 
EACH SUM 

PERIOD 

(HOURS) 

POWER (%) 
EACH SUM 

1 

12.68 

7.73 

7.73 

22.89 

10.47 

10.47 

2 

11.10 

5.76 

13.50 

19.13 

6.88 

17.36 

3 

11.55 

5.23 

18.73 

20.61 

6.57 

23.94 

4 

10.96 

5.08 

23.82 

13.68 

4.50 

28.44 

5 

10.26 

3.17 

27.00 

17.59 

4.43 

32.88 

6 

11.88 

2.97 

29.98 

16.21 

4.41 

37.30 

7 

8.79 

2.93 

32.91 

12.85 

4.12 

41.43 

8 

9.67 

2.61 

35.53 

68.14 

2.53 

43.97 

9 

10.37 

2.54 

38.07 

10.96 

2.47 

46.44 

10 

10.13 

2.26 

40.33 

12.05 

2.42 

48.87 

11 

9.26 

2.23 

42.57 

28.20 

2.09 

50.96 

12 

10.52 

2.16 

44.74 

29.20 

2.06 

53.03 

13 

14.79 

2.10 

46.84 

11.55 

2.01 

55.05 

14 

9.49 

1.99 

48.84 

10.52 

1.98 

57.03 

15 

12.05 

1.94 

50.79 

31.16 

1.81 

58.84 

16 

12.85 

1.90 

52.69 

9.49 

1.79 

60.63 

17 

13.68 

1.88 

54.57 

15.51 

1.42 

62.06 

18 

8.42 

1.74 

56.32 

12.68 

1.41 

63.48 

19 

20.61 

1.70 

58.02 

9.26 

1.39 

64.87 

20 

10.02 

1.65 

59.68 

43.99 

1.27 

66.15 


The three most powerful modes in the M2 spectrum are the modes with periods of 12.68 hours (7.73%), 11.10 hours 
(5.76%) and 1 1.55 hours (5.23%). For the K1 spectrum the three most powerful modes have periods of 22.89 hours (10.47%), 
19.13 hours (6.88%) and 20.61 hours (6.57%). Figures 1 to 6 below give the structure (amplitude and phase) for these 6 
modes. The amplitudes have been normalized to a maximum value of 100; the amplitude contours are shown by the solid 
lines and the dotted lines correspond to the phases. 

Figures 7 and 8 show the M2 and K1 tidal solutions. These solutions were obtained by including the first 150 eigen- 
functions from the <p and \jj solutions. The forced solution for 77 then contains 150 coefficients R^. The contours of equal 
amplitude and phase are given by the solid and dotted lines respectively. The arrows indicate the sense of progression for 
high and low water. The M2 solution has an amplitude RMS of 53.4 cm; the corresponding value for K1 is 10.0 cm. 
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OF POOR QUALITY 


























Fitting the M2 Solutions by Mazzega and Schwiderski 

The first effort to estimate the deep ocean tides from satellite altimeter data was done by Mazzega (1985). He used the 
SEASAT data to obtain a spherical harmonic expansion for the M2 tide in the world oceans. Mazzega’s solution reproduces 
many of the features predicted by the present numerical models, therefore a fit to his solution using Proudman functions 
should be of some interest. We have used Mazzega’s spherical harmonics coefficients to create values of amplitude and phase 
for the M2 tide at the 1915 points of our grid. We have used his expansion to degree and order 6; the resulting tide is shown 
in Figure 9 below. 



Figure 9. Mazzega’s M2 Tidal Solution 


The RMS in sea surface height -difference between our theoretical solution (150 functions) and Mazzega’s solution is 39.06 
cm before the fit. We then used the 1915 data points to estimate 150 coefficients which were used to create a fit to Mazzega’s 
solution. The RMS in sea surface height-difference between the fit and Mazzega’s solutionis 3.91 cm. Letting 1915 represent 
100%, the amplitude results are 539 points (28.14%) with fits within 5%, 1054 (55.03%) with fits within 10% and 1525 
points (79.63%) with fits within 20% of the Mazzega’s amplitude. The results for the phase are 650 points (33.94%) within 
5 degrees, 1425 points (74.41%) within 10 degrees and 1749 points (91.33%) within 20 degrees of the phase values given by 
Mazzega. Figure 10 shows the contours of amplitude difference (cm) between the fit and Mazzega’s solution. 

Schwiderski’s M2 solution was obtained from a tape where the values of amplitude and phase are specified on a 1 x 1 
grid; this grid is not congruent with the 3° x 3° grid used in the generation of the Proudman functions; there is a separation 
of a quarter of a degree in latitude and longitude between the points in the two grids. This separation was considered small 
for the purpose of our test and the values of amplitude and phase at the closest point of the 1° x 1° grid were simply 
adopted as the values in the 3° x 3° grid. Figure 1 1 below shows the result. 
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Figure 10. Amplitude Difference After Fitting Mazzega’s M2 Tide with 
150 Proudman’s Functions 











The RMS in sea surface height-difference between our theoretical solution and Schwiderski’s solution is 33.93 cm before 
the fit. A fit to Schwiderski’s solution is then obtained by estimating 150 coefficients from the 1915 data points. The RMS 
in sea surface height-difference between the fit and Schwiderski’s solution is 6.53 cm. The amplitude results are as follows: 
869 points (45.3%) with fits within 5%, 1356 points (70.8%) with fits within 10% and 171 1 points (89.3%) with fits within 
20% of Schwiderski’s amplitude. The phase results are 1479 points (77.2%) within 5 degrees, 1705 points (89.0%) within 10 
degrees and 1826 points (95.3%) within 20 degrees of the phase values given by Schwiderski. Figure 12 shows the contours 
of amplitude difference (cm) between the fit and Schwiderski’s solution. 

Spherical harmonics were also used to compute fits to Schwiderski’s M2 solution. For an expansion of degree = 7, invol- 
ving 128 coefficients, the RMS in sea surface height difference is 10.7 cm. The fit gets better with increasing degree, for an 
expansion of degree - 14 (450 coefficients) the RMS is down to 7.36 cm. Figure 13 shows the contours of amplitude dif- 
ference (cm) between this last case and Schwiderski’s solution. 

An examination of Figures 12 and 13 indicates that most of the error for both fits occurs in the boundaries of the ocean 
basin, certain areas being more critical than others. The highest errors occur in the Labrador Sea, the Northeast Atlantic, 
the Mozambique Channel, the Bay of Bengal, the Timor Sea, the Weddell Sea and in the spherical harmonic fit: the Patagonian 
Shelf. The artificial closure along the Tasman Plateau shows little error in both cases. 



Figure 12. Amplitude Difference After Fitting Schwiderski’s M2 Tide 
with 150 Proudman’s function 
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Figure 13. Amplitude Difference After Fitting Schwiderski’s M2 Tide with 
a 14th-degree Spherical Harmonic Expansion 


5. CONCLUSIONS 

The stream function and velocity potential orthogonal functions were computed numerically in the Atlantic-Indian 
Ocean basin using a 3° x 3° finite difference grid and a Lanczos procedure. These solutions were used to obtain the normal 
modes of the basin and the results were compared with those obtained previously using a 6° x 6° grid. The Laplace Tidal 
Equations were used to obtain the forced solutions for the M2 and K1 tidal components; these solutions do not include 
friction, self-attraction of the tide or tidal loading. A power spectrum for the M2 and K1 components was obtained, show- 
ing which modes are the most energetic. 

The eigenfunctions of the velocity potential were used to obtain fits to both Mazzega’s and Schwiderski’s solutions for 
the M2 tidal component. Using 150 functions Mazzega’s solution was fit with RMS = 3.91 cm in sea surface height-differ- 
ence. Schwiderski’s solution was fit with RMS = 6.53 cm; this is still lower than a 14th degree fit using spherical harmonics 
which yields RMS = 7.36 cm and involves three times the number of coefficients (450). 
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